# Dickstein, Ho, and Mark (2023)
# This script creates a output statistics related to counterfactuals results
# as we loop through different sets of parameters.

# Preliminaries
setwd("../../../../library")
source("PreliminariesCode.R")

## Which spec do you want to use?
spectouse <- "lowrisk_simplemh_censor0.2"

## Where is the counterfactual folder?
counterfactual_folder <- paste0(project_folder, "/analysis/counterfactuals")

# where is datasets folder
data_folder <- paste0(project_folder, "/data")

# What alpha-psi cutoff to use: 
alphampsicutoff <- .05

# What are the parameter sequences?
lambdashifter_sg <- rep(seq(.73, 1, .03), times = 10, each = 1)
lambdashifter_ind <- rep(seq(.73, 1, .03), times = 1, each = 10)
alphashifter_sg <- rep(1, 100)
alphashifter_ind <- rep(1, 100)
for(i in 1:length(lambdashifter_sg)){
  alphashifter_sg[i] = 1/lambdashifter_sg[i]
  alphashifter_ind[i] = 1/lambdashifter_ind[i]
}

# Load data
load(paste0(project_folder, "/analysis/counterfactuals/loop_specs/output/output_5_results_test.Rdata"))
load(paste0(project_folder, "/analysis/counterfactuals/loop_specs/output/output_base_results_test.Rdata"))
load(paste0(project_folder, "/analysis/counterfactuals/loop_specs/output/output_param_mat_test.Rdata"))

# Create table of end parameters
table_list <- list()
for(i in seq_along(results_list_5)){
  table_list[[i]] <- data.table(
    im.beta1.sum_concurrent_risk = param_mat$im.beta1.sum_concurrent_risk[i], 
    im.beta2 = param_mat$im.beta2.cons[i], 
    im.beta3 = param_mat$im.beta3.cons[i],
    sg.beta1.sum_concurrent_risk = param_mat$sg.beta1.sum_concurrent_risk[i], 
    sg.beta2 = param_mat$sg.beta2.cons[i], 
    sg.beta3 = param_mat$sg.beta3.cons[i],
    markup = param_mat$mrkup[i],
    im.post.cs_1 = results_list_5[[i]][[1]]$welfare["CS_1"]*100,
    im.post.cs_2 = results_list_5[[i]][[1]]$welfare["CS_2"]*100,
    im.post.g_spend = results_list_5[[i]][[1]]$welfare["GovmntSpend"], 
    im.post.unin = results_list_5[[i]][[1]]$mktshr["PUninsured"], 
    im.post.bron = results_list_5[[i]][[1]]$mktshr["PBronze"], 
    im.post.silv = results_list_5[[i]][[1]]$mktshr["PSilver"], 
    im.post.gold = results_list_5[[i]][[1]]$mktshr["PGold"], 
    im.post.gross_p_b = results_list_5[[i]][[1]]$grossprem[1], 
    im.post.gross_p_s = results_list_5[[i]][[1]]$grossprem[2],
    im.post.gross_p_g = results_list_5[[i]][[1]]$grossprem[3],
    im.post.net_p_b = results_list_5[[i]][[1]]$netprem[1],
    im.post.net_p_s = results_list_5[[i]][[1]]$netprem[2],
    im.post.net_p_g = results_list_5[[i]][[1]]$netprem[3]*100,
    im.pre.cs_1 = results_list_base[[i]][[1]]$welfare["CS_1"]*100,
    im.pre.cs_2 = results_list_base[[i]][[1]]$welfare["CS_2"]*100,
    im.pre.g_spend = results_list_base[[i]][[1]]$welfare["GovmntSpend"], 
    im.pre.unin = results_list_base[[i]][[1]]$mktshr["PUninsured"], 
    im.pre.bron = results_list_base[[i]][[1]]$mktshr["PBronze"], 
    im.pre.silv = results_list_base[[i]][[1]]$mktshr["PSilver"], 
    im.pre.gold = results_list_base[[i]][[1]]$mktshr["PGold"], 
    im.pre.gross_p_b = results_list_base[[i]][[1]]$grossprem[1], 
    im.pre.gross_p_s = results_list_base[[i]][[1]]$grossprem[2],
    im.pre.gross_p_g = results_list_base[[i]][[1]]$grossprem[3],
    im.pre.net_p_b = results_list_base[[i]][[1]]$netprem[1],
    im.pre.net_p_s = results_list_base[[i]][[1]]$netprem[2],
    im.pre.net_p_g = results_list_base[[i]][[1]]$netprem[3]*100,
    sg.post.cs_1 = results_list_5[[i]][[2]]$welfare["CS_1"]*100,
    sg.post.cs_2 = results_list_5[[i]][[2]]$welfare["CS_2"]*100,
    sg.post.g_spend = results_list_5[[i]][[2]]$welfare["GovmntSpend"], 
    sg.post.unin = results_list_5[[i]][[2]]$mktshr["PUninsured"], 
    sg.post.bron = results_list_5[[i]][[2]]$mktshr["PBronze"], 
    sg.post.silv = results_list_5[[i]][[2]]$mktshr["PSilver"], 
    sg.post.gold = results_list_5[[i]][[2]]$mktshr["PGold"], 
    sg.post.gross_p_b = results_list_5[[i]][[2]]$grossprem[1], 
    sg.post.gross_p_s = results_list_5[[i]][[2]]$grossprem[2],
    sg.post.gross_p_g = results_list_5[[i]][[2]]$grossprem[3],
    sg.post.net_p_b = results_list_5[[i]][[2]]$netprem[1],
    sg.post.net_p_s = results_list_5[[i]][[2]]$netprem[2],
    sg.post.net_p_g = results_list_5[[i]][[2]]$netprem[3]*100,
    sg.pre.cs_1 = results_list_base[[i]][[2]]$welfare["CS_1"]*100,
    sg.pre.cs_2 = results_list_base[[i]][[2]]$welfare["CS_2"]*100,
    sg.pre.g_spend = results_list_base[[i]][[2]]$welfare["GovmntSpend"], 
    sg.pre.unin = results_list_base[[i]][[2]]$mktshr["PUninsured"], 
    sg.pre.bron = results_list_base[[i]][[2]]$mktshr["PBronze"], 
    sg.pre.silv = results_list_base[[i]][[2]]$mktshr["PSilver"], 
    sg.pre.gold = results_list_base[[i]][[2]]$mktshr["PGold"], 
    sg.pre.gross_p_b = results_list_base[[i]][[2]]$grossprem[1], 
    sg.pre.gross_p_s = results_list_base[[i]][[2]]$grossprem[2],
    sg.pre.gross_p_g = results_list_base[[i]][[2]]$grossprem[3],
    sg.pre.net_p_b = results_list_base[[i]][[2]]$netprem[1],
    sg.pre.net_p_s = results_list_base[[i]][[2]]$netprem[2],
    sg.pre.net_p_g = results_list_base[[i]][[2]]$netprem[3]*100,
    sg.lambdashifter_sg = ifelse(i >= 234, rep(seq(.73, 1, .03), times = 10, each = 1)[i-233], 1),
    im.lambdashifter_ind = ifelse(i >= 234, rep(seq(.73, 1, .03), times = 1, each = 10)[i-233], 1)
  )
}
table <- do.call("rbind", table_list)

psi_upperbound <- .3

# Adding in small group consumer surplus pre 
sg_surplus <- fread(paste0(counterfactual_folder, "/specs/", spectouse, "/data/SGSurplusEstimationFile.csv"))
## dropping outside options and keeping 2016
sg_surplus <- sg_surplus[!(x_av == 0)]
sg_surplus <- sg_surplus[year == 2016]
sg_surplus <- sg_surplus[!(fixedeffect_beta0 == 0 & x_av != 0)]

# Generating surplus:
print("generating small group pre surplus:")
for(i in seq_along(results_list_5)){
  sg_surplus[, omega := exp(param_mat$sg.beta2.cons[i]), ]
  sg_surplus[, psi := exp(param_mat$sg.beta3.cons[i]), ]
  sg_surplus2 <- sg_surplus
  
  if(i >= 234 & i <= 333){
    sg_surplus2[, in_orig_cs1 := ifelse(alpha - psi > 0, 1, 0)]
    sg_surplus2[, in_orig_cs2 := ifelse(alpha - psi > alphampsicutoff, 1, 0)]
    sg_surplus2$alpha <- sg_surplus$alpha*alphashifter_sg[i-233]
  } 
  
  sg_surplus2[, V_notiering := x_av + .5 * omega * (x_av ^ 2) - (alpha - psi) * price_hh_notiering + fixedeffect_beta0, ]
  sg_surplus2[, V := x_av + .5 * omega * (x_av ^ 2) - (alpha - psi) * price_hh + fixedeffect_beta0, ]
  
  if(i > 81 & i <= 181){
    sg_sur_collapsed <- sg_surplus2[x_av > 0, .(
      logsum = log(sum(exp(V))), 
      alphampsi = mode1(alpha - psi_upperbound)),
      by = c("subscriberid", "year")]
  } else if(i >= 234 & i <= 333) {
    sg_sur_collapsed <- sg_surplus2[x_av > 0, .(
      logsum = log(sum(exp(V))), 
      alphampsi = mode1(alpha - psi), 
      in_orig1 = sum(in_orig_cs1), 
      in_orig2 = sum(in_orig_cs2)),
      by = c("subscriberid", "year")]
    
  } else {
    sg_sur_collapsed <- sg_surplus2[x_av > 0, .(
      logsum = log(sum(exp(V))), 
      alphampsi = mode1(alpha - psi)), 
      by = c("subscriberid", "year")]
  }
  
  sg_sur_collapsed$cs <- sg_sur_collapsed$logsum / sg_sur_collapsed$alphampsi
  if(i >= 234 & i <= 333){
    table$sg.pre.cs_1[i] <- mean(sg_sur_collapsed[(alphampsi > 0) & (in_orig1 >= 1)]$cs)*100
    table$sg.pre.cs_2[i] <- mean(sg_sur_collapsed[(alphampsi > alphampsicutoff) & (in_orig2 >= 1)]$cs)*100
  } else {
    table$sg.pre.cs_1[i] <- mean(sg_sur_collapsed[alphampsi > 0]$cs)*100
    table$sg.pre.cs_2[i] <- mean(sg_sur_collapsed[alphampsi > alphampsicutoff]$cs)*100
  }
  print(paste0(i))
}

# Saving
fwrite(table, paste0(data_folder, "/outcome_table_hm.csv"))
